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Abstract. The non-linear evolution of a stratified perturbation in a three dimensional expanding Universe is 
considered. A general Lagrangian scheme (Q model) is introduced and numerical investigations are performed. 
The asymptotic contraction of the core of the agglomeration is studied. A power-law scaling is detected and 
an heuristic interpretation of the numerical findings is provided. An asymptotic equation for the multi-stream 
velocity flow is derived and it is shown to agree quantitatively with the dynamics of the Q model. The relation 
to the adhesion model is discussed. 
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1. Introduction 

The present Universe is inhomogeneous with structures 
of many scales, from galaxies to galaxy clusters and su- 
perclusters. Around 100000 years after Big Bang, the 
Universe was however very nearly homogeneous, with den- 
sity fluctuations of relative magnitude of about 1 : 10 5 . 
The presently observed large scale structure have been 
generated by the process of the gravitational instability, 
acting on those initially small perturbations. 

The goal of this paper is to investigate in direct numer- 
ical simulations and theoretical analysis the validity of the 
adhesion model, a benchmark for the non linear develop- 
ment of the gravitational instability. We will make con- 
tact with recent results of Buchert & Dominguez ( 1999 ) 
and present a sharpening of their main findings. The con- 
clusion of this study is that the Burgers equation of the 
adhesion model does not model structure formation quan- 
titatively correctly, but another transport equation of a 
similar type does. 

We will limit ourselves to a flat, critical Universe 
(fi = 1) where most mass is contained in a dark com- 
ponent. In fact, we assume effectively that all the matter 
is dark, or behaves as such. There is no vacuum energy, 
or, equivalently, the cosmological constant is zero. This 
has been the favored cosmological model in the recent 
past and, although it is not at the present, it is still of 
conceptual and qualitative interest. We remark that the 



approximation of a flat critical Universe will be an accu- 
rate one for most of the time since recombination epoch, 
also in the presently favored models. 

We will further make the assumption of stratified per- 
turbations. This is surely not realistic, but, as we will see, 
it allows for a numerical scheme of incomparable speed 
and accuracy, describing the full process of structure for- 
mation, and is therefore a useful testing ground. 

The paper is organized as follows: in Section ^| we 
present the general background. In Section ^ we recall the 
derivation of the Quintic model (Aurell & Fanelli |200ip , 
that allows to study the evolution of a stratified perturba- 
tion in a Einstein-de Sitter, three dimensional expanding 
Universe. Section [| is devoted to the discussion of the nu- 
merical implementation. In Section^ we present the result 
of our numerical study: scaling laws are displayed and an 
heuristic explanation is provided. In Section |fij we derive 
an equation of transport that is consistent with the nu- 
merical findings. In the final Section |^ we sum up and 
discuss our results. 



2. General Background 

The linear regime of structure formation was studied by 
Lifshitz 1947 and is described in several classical mono- 
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graphs (Peebles 1980, Weinberg 1972). For the non lin- 
ear regime, there is a long tradition of considering various 
simplified models, from the Zeldovich pancake (Shandarin 
& Zeldovich 1989) and the adhesion model ( Gurbatov, 
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Saichev & Shandarin 1989 ), to full-scale iV-body simula- 
tions, with the accompanying approximations of numeri- 
cal nature. We start with a short review of the basic setup 
and equations. 

Collision-less dark matter is described by the kinetic 
Vlasov-Poisson equations, in an expanding three dimen- 
sional Universe. This level of approximations contains 
the assumption that tensor degrees (gravitational waves) 
can be neglected, and that the matter motion is quasi- 
Newtonian. We assume an inertial reference frame, and 
label the position with r. Then it is customary to intro- 
duce the comoving coordinate x by performing the follow- 



ing transformation (Peebles 1980) 



a(t)x , 



(1) 



where the scale factor a(t) is function of proper world time. 
For the critical Friedman Universe: 



" = I — 



2/3 



(2) 



and t 2 — 67rGp(io), where p(io) is the homogeneous den- 



sity at time t (Peebles |1980| , Weinberg |1972| ) and G is 
the gravitational constant. The Vlasov-Poisson equations 
read: 



d t f + 



V x / - V x ^ • V p / = 



V 2 4> = 4nGa 2 (p - 



Pb) 



where p is the variable conjugated to x; /(x, p, t) is the 
distribution function in the six dimensional phase space 
(x, p) ; ifi is the gravitational potential and pb is the mean 
mass density distribution. The particle density p(x, t) and 
velocities u(x, t) are given, in term of /(x, p, t), as: 



TTi f 

p(x,t) = ^3 / /(x,p,i)dp , 

p(x, i)u(x, t) = ^ J p/(x,p,t)dp. 



(4) 



(5) 



It is well known that (|^) admits special sol utions of 
the form (Vergassola, Dubrulle, Frisch & Noullez 1993 ): 



/(x,p,i) = 



S (p — mau(x, t)) , 



(6) 



where d is the dimension of space and S d (.) the d- 
dimcnsional delta function. We will refer to this class as 
single-speed solutions, because to each given (x, t) corre- 
sponds a well defined velocity u. Assuming (||), after some 
manipulation, it follows from equations (Q) and (||): 



d t p + 3-p+-V ■ (pu) = 
a a 



_ a 1 . . 

o t u H — u H — (u • Vju = g 

a a 



, V • g = -4irGa(p - p b ) 



(7) 



where we have introduced g = ~Vip/a such that V x g = 
0. It should be stressed that system (Q), is valid as long 
as the distribution function /(x, p) is in the form (^J). 
Beyond the time of caustic formation, when the fast parti- 
cles cross the slow ones, the solution become multi-stream. 
Hence, the pressure-less and dissipation-less hydrodynam- 
ical equations (uh are incomplete. We will in this paper 
extend system (R) beyond caustic formation for strati- 
fied flows, i.e. when velocity has one component only, and 
varies with respect to this direction. 

Let us focus now on time before the first particle cross- 
ing. Then we can further make the assumption of par- 
allelism: the peculiar velocity is a potential field, which 
remains parallel to the gravitational peculiar accelera- 
tion field (B uchert T., Dominguez A. & Perez-Mercader 
1999 L Peebles [l980l Vergassola, Dubrulle, Frisch & Noullez 
1993|) : 

(8) 



g = F(t)u 



where F(t) is a positive, time dependent, proportional- 
ity coefficient. Relation (||) is well justified in the linear, 
as well in the weakly non linear regimes and allows to 
treat analytically the problem. From the linear theory 
it fol lows (B uchert T., Dominguez A. & Perez-Mercader 
1999L Peebles [l980l Vergassola, Dubrulle, Frisch & Noullez 
1993|) : 



(3) F(t) = 4irGp b b/b 



(9) 



where b represents the growing mode of the density field 
in the linear regime. Hence, defining the new velocity field 
v = u/(ab), the system M) reduces to: 



d b v + (v • V)v = 



(10) 



which is the multidimensional Burgers equation 
(Vergassola M., Dubrulle B., Frisch U. & Noullez A. |1993[ ). 
The inviscid Burgers equation describes the free motion 
of fluid particles subject to zero forcing and is equivalent 
to the famous Zeldovich approximation (Shandarin & 
Zeldovich 1989). Again, it is worth recalling that the 



picture is correct as long as the solution stays single- 
stream. After caustic formation, it has been proposed 
that the resulting change on the gravitational force can 
be modeled by an effective diffusive term (adhesion model 
(Gurbatov S.N., Saichev A.I. & Shandarin S.F. |l989| )). 
This should represent the effect of the gravitational 
sticking not captured by the Zeldovich approximation. 
Mathematically, this means introducing a term of the 
form /xV 2 v, in the right hand side of the equation (|l0|). 
In order for the diffusion term to have a smoothening 
effect only in those regions where the particles crossing 
takes place, the viscosity p, should be small. The adhesion 
model reads: 



d b v + (v • V)v = pV 2 v 



< v = -VV> 



k d b p + V • (pv) = 



(11) 
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where ip — tp/(bF(t)). The limit when viscosity p tends to 
zero is often taken. We remark, that, as is well known this 
is not equivalent to setting p to zero from the outset, but 
is instead a regularization of (|l0|), equivalent to the Lax 
entropy condition. 

Although numerical experiments suggest qualitative 
agreement, no theory is to our knowledge presently avail- 
able that quantifies the exact relationship between (^) and 
(|li"|), after caustic formation. 

The hypothesis has been put forward that (^lj) is an 
asymptotic description of (||) , for particular classes of ini- 
tial conditions (Starobinsky 2000). This statement relies 
on the assumption that, because of the expansion, the 
thickness of the individual pancakes grows more slowly 
than their typical separation. We will refer to this picture 
as to Starobinsky conjecture. Repeating the statement in 
the introduction, the aim of this paper is to investigate 
the asymptotic relation between the Vlasov approach (|J) 
and the adhesion model ([ll]), focusing our attention on 
the stratified dynamics. To jump also again to the conclu- 
sion, we will derive transport equation similar to those of 
Buchcrt& Dominguez (1999), but not quite the same. 



3. Quintic model 

We will study stratified perturbation in a critical Universe, 
in a particle representation. We recall that a non-linear 
change of time allows us to transform the equations of 
motion to ordinary differential equations with constant 
coefficients, and which can therefore be integrated by a 
fast event-driven scheme. We will refer to th is dynamics 
as to the Quintic model (Aurell & Fanelli 2001 ) , for reason 
that will become clear. 

The Newtonian equations of motion for N particles 
follow from a Lagrangian 



(12) 



where ~V 2 (f> = AirGp. In the point particle picture the den- 
sity profile is: 



(13) 



where Xi is the comoving coordinate of the i — th particles, 
in the direction of which the density and velocities varies. 

Expressing ([l2]) as function of the proper coordinate, 
Xi, and assuming (|l3|), the equation of motion of the i-th 
particle reads: 

d 2 x ' Qj dx ■ 

-^ + 2-— --^Gp b {t) Xl =a- 3 E grav { Xl ,t) , (14) 
at a dt 

where 



Egrav{xi,t) = -2-kG mjsign(xj - Xj) 

3 

From the equation of continuity: 
pb(t) = poa(t)~ 3 , 



(15) 



(16) 



and, by performing a suitable non linear transformation 
of the time variable it is possible to concentrate all the 
time dependence in the term 2^^-. The choice is: 



dt = a 3/2 dr , 



(17) 



where r has dimension of time ( Rouet, Feix & Navet 1990 



Rouet et al. 1991). The equation (14) is thus transformed 
into: 



d 2 Xi 



2 dr 



AirGpoXi = E grav (xi, t) . 



(18) 



In a flat Einstein de Sitter model the scale factor a(t) 
grows with time as a power- law (see eq. (^)) and therefore 
cq. (HI) takes the form: 



d 2 Xi 1 dxi 
dr 2 3io dr 



g ,2 X — E grav yX$ , 7") 
6t 



Q model . (19) 



The possibility of making this coordinate change is in- 
deed the technical reason why we work with a critical 
Universe. Equation fll9| ) is the model we refer to as to 
the Quintic (Q) model. We stress that the quintic model 
is nothing but a particle picture of the full self-gravitating 
dynamics for the class of stratified perturbations. The con- 
tinuum (N — ► oo) limit of ( fi9| ) is just (||), with different 
time coordinate. The interest of this formulation is that, 
as for the classical static self-gravitating systems in one 
dimension, E grav is a Lagrangian invariant, proportional 
to the net mass difference to the right and to the left of 
a given particle at a given time. Hence, the evolution of 
the system is recovered by identifying the time and loca- 
tion of the particles crossings, and connecting analytical 
solut ions between such events (Noullez, Fanelli & Aurell 

teooil) . 



4. Numerical scheme 

The equation of motion of each particle, in the Q model, 
is specified by equation (|l^). In between successive cross- 
ings (t 6 [t",t" +1 ]) the right hand side is a constant 
and, therefore, ( |l9|) adm its an explicit solution in the form 
(Aurell & Fanelli |200l|): 



Xi(r) = c\ exp( 



2(7 



3t 



') (r 
-) + c\ cxp(-^- 



to 



where K t = (3tl/2)E grav (xi,T), is constant. The coef- 
ficients c\ and c\ are determined by xf = Xi(r n ) and 
wf = Xi(r n ), i.e. by the states of the particle at the time 
of the last crossing, and read: 



tow? 



2 r 



Ki 



Ki 



(21) 



The form of equation (j20j) suggests introducing an aux- 
iliary variable z = exp((r — r n )/3to). The crossing times 
between neighboring particles (i.e. i,i + 1) can, hence, be 
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computed by solving numerically a quintic equation in the in Fig. |l], the only real root of eq. (22) lies in the interval 



form: 

/(*) = Al l+1 z 5 - B^ +1 z 3 + C? M1 = , 
where: 



(22) 



All 



BV 



-[Axf+toAvf-iKw-Ki)] 



-(Ki 



Ki 



const 



(23) 



°M+l - 5 



Ax" - -t Aw? 



K,. 



and Ax™ = xf +1 — xf, Aw" = wf +1 — w™. The function 
f(z) represents the distance between particles i and i + 1, 
at transformed time z. 

The event-driven scheme (Noullez, Fanelli & Aurell 



2001) can be adopted to follow the dynamics of the system. 
The positions of the particles are stored in monotonically 
increasing order. A proper time, 6i, is associated to each 
particle i: it refers back to the time the particle last ex- 
perienced a collision. Initially all Qi are set to zero. The 
algorithm computes first the crossing time of each particle 
with the neighbor to the right, by solving N — 1 quintic 
equations, as in ( |22| ) (see above). The results are then 
stored in an array, which is sorted on a heap. Once the 
heap has been built, the minimum collision time, t m i n , is 
at the root, i.e. at the first position in the heap. The parti- 
cles involved in the first collision are selected by means of a 
trivial 0(1) operation. The algorithm lets them evolve up 
to t m i n , according to the equation (|2C|). At this point, the 
particles share the same spatial position, and the crossing 
takes place (the velocities are exchanged). The successive 
step is to compute the next predicted collision time be- 
tween i and i + 1. The new value replaces the old one, 
and the heap needs to be rearranged. In addition, as an 
effect of the changes in velocity for particles i and i + 1 
particles, the two collisions with their nearest neighbors 
(respectively i — 1 and i + 2) need to be re-computed. The 
heap is then re-arranged with at most (log (AT)) opera- 
tions and the whole procedure can be repeated for N co u 
collisions. Therefore the complexity of the algorithm is 
in wo rst-case 0(N co u log(AT)) (Noullez, Fanelli & Aurell 
200ll ). 

Indeed, in this particular application, the numerical 
procedure for finding the solution of the quintic equation 
might converge slowly, affecting the whole computation 
time, with a non negligible contribution. Therefore, in or- 
der to achieve the goal of a fast implementation, special 
care has to be devoted to the analysis of (p2|). Recalling 
the definition of f(z), we look for the smallest real value 
z > 1, such that f(z) = 0. As a first trivial observation, 
we notice that, by definition, z is larger than one, and the 
inter-particle distance f(z = 1) is non- negative. In addi- 
tion, B™ i+1 is a positive constant, independent of i. 

Two scenarios are therefore possible: if the coefficient 
Af i+i is positive, no crossing is allowed. As it is sketched 



[0,1], and therefore has to be rejected. This means that 




Fig. 1. The function f(z) represented versus z, for Af i+1 
positive. The solid line refers to the case when C™ i+1 is 
positive while the dashed line to C" i+1 negative. Note that 
no intersection with the horizontal axis is allowed for z > 
1 (since /(l) > 0, by definition). 



these two particles will never collide if left to themselves: 
the expansion is too strong to be overcome by their mutual 
gravitational attraction. 

On the other hand, if the coefficient A™ i+1 is nega- 
tive, more care is required. The problem is to bound the 
root in a reasonable interval in order to assure a fast con- 
vergence of a numeric procedure. First, we observe that 
there is a local maximum for z > 0. The coordinate z max 
is easily computed and f(z max ) is evaluated; f(z max ) is 
positive, since it is by definition larger of /(l) > 0. Since 
the function f(z) goes to — oo as z — > oo, there should be 
an intersection [z > 1 s.t. f{z) = 0) with the horizontal 
axis, see Fig ||. The following procedure is adopted. First 
we introduce h(z) (thin solid line in Figs. ||, ||), that is the 
quadratic approximation of f{z) around z max , defined by: 

^(^) — f (^mai) "T" ~/ (z max ^(z Z max ^ . (^^) 

Then we compute the intersection z\ of h(z) with the 
horizontal axis: 



zi 



f^max 

w 



-5 a: 



M+i 



3B 



(25) 



Two different cases have to be considered (respectively, 
Figs. ||||). If / (zi) is positive, we take the tangent to f(z) 
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Fig. 2. The function f(z) represented versus z, for A™ i+1 
negative (thick solid line). The thin solid line represents 
the quadratic approximation around z — z rnax . Here f(z\) 
is negative and, therefore, the root z is bounded between 

[Zmax? Z\\ • 




Fig. 3. The function f{z) represented versus z, for Af i+1 
negative (thick solid line). The thin line is the quadratic 
approximation around z = z max . Here f(zi) is positive 
and therefore, the root z is bounded between [z\ , Z2] ■ Here 
Z2 is the intersection with the z axis of the tangent to f(z) 
in z\ (dotted line). 



in Z\, and compute its intersection, z% , with the z axis 
(Fig. |): 



By construction z^> z and, therefore, z € \z\,zi\ If, 
on the contrary, f{z\) is negative, then z € [z max , Z\], see 
Pig. I 

In both cases we have confined the root in a narrow 
interval: there, f(z) is a monotonic decreasing function. 
Therefore, we can apply a combination of bisection and 



Newton-Rapson method (Press 1992). This hybrid algo- 
rithm assures a stable and fast convergence to the solution. 
In the present implementation we assumed a tolerance er- 
ror of 10~ 13 . 

5. Asymptotic scaling laws: numerical results and 
heuristic interpretation 

We simulate the dynamics of the Q model, by using the 
numerical scheme discussed above. We consider a sys- 
tem of N particles of mass m = 1/N, confined in a box 
of size L. Reflecting boundary conditions are assumed, 
which is equivalent to consider periodic perturba tion o f 
size 2L, with reflexion symmetry (Aurell & Fanelli 2001 ). 
We choose units such that 4irG is equal to one. Time is 
measured by the a-dimensional quantity (t/to). The unit 
of length is the spatial interval in which the particles are 
initially distributed (i.e. L — 1), and thus the initial den- 
sity po is set to one. 

In particular, we are interested in the late time evolu- 
tion of the system to better understand the validity of the 
Starobinsky conjecture, the aim of this analysis being to 
provide a quantitative test of the reliability of the adhe- 
sion model, as an effective description of the mechanism 
of large scale structure formation in the Universe. 

With this in mind, we consider the evolution of an indi- 
vidual cluster and investigate the process of gravitational 
collapse. We measure the progressive contraction of the 
inner region of the agglomeration, compared to the over- 
all expansion. This effect can be computed by the ratio 
between the width of region, Ax, that contains half of the 
whole mass of the system, centered around the position of 
maximum density, and the size of the perturbation L []. 
This quantity is then plotted, as function of the rescaled 
cosmological time t/t . 

Simulations are performed for two different classes of 
initial conditions. In both cases, we make the non restric- 
tive choice of considering a perturbation centered around 
our center of reference. 

First, the initial velocity is assumed to be a smooth 
function of position. As clearly expected, the system de- 
velops spiral in the phase space like the one displayed in 
Fig. ||. In Fig. || we plot Ax/L, versus t/to- The exper- 
iments are performed for different values of N. A clear 
power-law behavior: 



Ax 



(27) 



Z2 = Zi 



f(zi) 



5A™ i+1 z 1 + 3Bf i+1 z 1 



(26) 



is always displayed, regardless of the number of particles 
simulated. This implies that the result is not affected by 

1 In comoving coordinate L is a constant, which we have set 
to one. 
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Fig. 4. Velocity field versus positions, starting from a sin- 
gle speed initial condition (sinus wave). Here N — 4096 
andt/to — 6.1249064 x 10 4 . Reflecting boundaries are as- 
sumed. Positions and velocities are in arbitrary units. 



the discreteness of the representation, and can be assumed 
to hold in the continuum limit. The scaling is present over 
several decades and the best numerical fit gives the value 
a = —0.22. Then, a source of noise is introduced in the 



gle realization. They show complete agreement with (JsJ) 
. ft is worth stressing that these results are not sensitive 



10" 




10 



10" 



10" 

///„ 



10 10 



Fig. 6. Ax/L is plotted versus t/to- Here wq represents 
the amplitude of the initial smooth sinus to which the white 
noise, s wn , is superposed (i.e., w = WQsin(x) + s wn , see 
small inset). Different symbols refer to different values of 
N (see the legend). The scaling is consistent with the one 
derived from Fig. (here represented by the thick solid 
line). 




t/t n 



Fig. 5. Ax/L is plotted versus t/to- Different symbols re- 
fer to different values of N (see the legend). A clear power- 
law behavior is displayed. The best numerical fit gives ex- 
ponent a — —0.22. The small inset represent the initial 
condition in the phase space (x,w). 



initial condition: as shown in the small inset of Fig. ^| a 
white noise signal is generated and superposed to a sinus 
wave of amplitude wq. The thickness of the dense region 
is studied, following the line of the preceding discussion. 
Results are reported in the main plot of Fig. ||, for a sin- 



to the choice of measuring the interval that contains N/2 
particles. Any other finite fraction leads to the same con- 
clusions. 

In order to provide a full picture, we performed sim- 
ilar analysis for the velocities distribution. Consider Aw, 
such that N/2 particles have velocities in the interval 
[Aw/2, —Aw/2]. The ratio Aw/wo is plotted vs. t/to, 
where wo represents the amplitude of the initial smooth 
wave (see captions of Figs. [7 and [8]). Again, and for both 
the initial conditions considered here, a power-law behav- 



ior: 

Aw 

w 



(28) 



is clearly displayed (Figs. 0, |[). Here the best numerical 
fit gives (3 = -0.11. 

Assuming the occurrence of power-law behaviors, there 
is a simple, heuristic, explanation for deriving the correct 
value of the exponents, in agreement with the numeri- 
cal findings. As previously stated, the Quintic model is 
equivalent to the Vlasov-Poisson set o f equations in the 
continuous limit (Aurell & Fanelli 2001). Such a system is 
conservative and therefore the volumes in the phase space 
are conserved (Liouvillc theorem). Hence: 



AxAp = const , 



(29) 



where Ax and Ap represent, respectively, any space and 
momentum intervals in the comoving reference frame. 
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Fig. 7. Aw/w is plotted versus t/t . Here Wq represents 
the amplitude of the initial velocity perturbation (i.e., 
w = wq sm(x), see small inset). Different symbols refer to 
different values of N (see the legend). A clear power-law 
behavior is displayed. The best numerical fit gives exponent 
(3 = -0.11 
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Fig. 8. Aw/ wo is plotted versus t/to- Here wq represents 
the amplitude of the initial smooth sinus to which the white 
noise, s wn , is superposed (i.e., w = wo sin(x) + s wn , see 
small inset). Different symbols refer to different values of 
N (see the legend). The scaling is consistent with the one 
derived from Fig. (here represented by the thick solid 
line). 



Recalling that w = dx/dr and that a(t) = (t/to) 2 / 3 , 
cq. ( p9[ ) is transformed into: 

( t\- 1/3 

Ax Aw = ( — J . (30) 

Assuming that Ax and Aw scales in time as power- 
laws, respectively with exponents a and (3 it follows that: 

1 



a + (3 = 



(31) 



A second relation is needed to close the system. Let 
us consider the Hamiltonian H = K + V which describes 
the Quintic model, except the friction term. K is a normal 
kinetic term, quadratic in the velocities, and V is made of 
two terms: the interaction and the background potential. 
The following relation holds: 



dt 6to 



(32) 



Assume that the global virialization has occurred. 
Then < K >~< V >, where < K > and < V > represent 
the time average of the kinetic and potential energies. K 
is quadratic in velocities while V has one term quadratic 
(the background), and one linear (the interaction). If we 
have a mass concentration at the origin, the quadratic can 
be ignored compared to the linear. 

Hence, neglecting the time averages, we have 
velocities 2 ~ distances. If Ax ~ (t/to) a , therefore Aw ~ 
(t/to) a / 2 , which implies: 

P = f , (33) 

and therefore a = —2/9 and [3 — —1/9, in agreement with 
the numerics. From eq. (p2[), the dissipative time is: 



Tdiss — H/(dH/dt) ~ constant 



(34) 



The virialization time of the mass agglomeration, T V i r 
is about the time it takes one particle to traverse the mass, 
that is 

( t V 

T v ir = (distance) / (velocity) ~ I — J . (35) 

Hence, the virialization time becomes quickly smaller 
than the dissipative time, and the argument is self- 
consistent. 

6. Transport equation 

As already observed, relation (||) holds as long as the so- 
lution stays single stream. Hence the hydrodynamical pic- 
ture (0) is applicable just before the time of caustic for- 
mation, when the first particles crossing takes place. In 
order to extend the analysis beyond that time we have to 
consider the more general solution to (w): 



/(x,p,t) 



a 3 p(x, t) 



/o(x,p,t) 



(36) 



where /b( x > P> t) is the velocity profile. 

Assuming ( p6[ ) and performing the same analysis as 
described in Section^, the following system is derived from 
Vlasov-Poisson equations (m : 



d t p + 3-p+-V-(pu) = 

a a 



< d t u+ -u+ - (u- V)u = g V 

a a ap 



p(u 2 



u 2 ) 



(37) 



V • g = -4irGa(p - p b ) 
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where [ • ] represents the mean, after averaging over ve- 
locity space. Comparing with (Q), we see that in the sec- 
ond equation of (|3^) an extra-term appears, that takes 
into account the effect of the dispersion in velocities, in 
the region where the multi-values solution is developed. 
When the flow is single-stream, the variance is zero and 
the system (37) is identical to (^) (u = u). This deriva- 



tion was already discussed in (Buchert & Dominguez 1999) 
and represents the starting point of our analysis. We will 
limit ourselves to the case of stratified perturbations, and 
derive an equation for the transport of the velocity flow, 
that agrees with the numerical analysis performed at the 
level of the Q model (i.e. discrete Vlasov). In this respect, 
we are mainly concerned by the progressive collapse of 
the inner bulk (see Figs a key feature that needs to 
be considered explicitly in our analysis. The question that 
naturally arises is whether or not the transport is well rep- 
resented by Burgers' equation. Here on, we will indicate 
the mean velocity simply with u. 

Let us consider the dynamics of the Q model. It can be 
shown numerically that, for all classes of initial conditions 
here considered, far beyond the time of first crossing, the 
multi-stream velocity profile fo(w) is well approximated 
by a Gaussian: 



fo{w) oc exp — ) , 

T(x) 



(38) 



where the temperature , T(x), is smaller in the bulk (more 
narrow Gaussian) than in the surrounding regions (larger 
profile). Equation (38) is also, of course, the natural choice 
in a system close to equilibrium, where p,w and T vary 
considerably over a cluster, but little over distances of the 
order of inter-particle spacing. 

Our ansatz for T(x) is the following: 



T(x) = K 



(t/t 



(39) 



where K is a constant and 7 is a real number belonging 
to the interval [0, 1]. This choice reproduces qualitatively 
the shrinking of the velocity profile, observed in comoving 
coordinates in correspondence of the denser inner regions. 
On a more quantitative level, numerical checks have been 
performed to support the validity of (|39|), showing in all 
cases, a satisfactory agreement. In particular, from ( pSj ) it 
follows that the variance, var^(x), is given by: 



var®(x) = uu 2 



w 



2 p-y 



(40) 



where the label Q indicates that we work in Quintic vari- 
ables. In the main plot of figure (^|) var^(x) is plotted ver- 
sus x. The filled circles refers to the results of the numer- 
ical experiments, while the solid line is a two parameters 
fitting, based on (p9|). (i.e. var®(x) — A\/p A2 , where now 
p is the histogram of positions and A\ , Ai free parame- 
ters). In Fig. Ai = 0.3: the same value, within the errors 
(~ 0.1), is found for different initial realizations and/or 
different times. Larger deviations from (E3) are observed 



CO 

> 

CD 
O 

c 

CO 
'1- 
CO 

> 



0.03 



0.02 



0.01 - 




-0.07 



-0.035 



0.035 



0.07 



positions, x 



Fig. 9. Main plot: var® (x) is plotted versus x. Here 
N = 8192 and t/t = 6.1249064 x 10 4 . The filled circles 
refers to the numerical experiment. The dashed line is a 
guide to the eye. The solid line is obtained by making use 
of the fitting function A\/ p A2 , with A\, Ai, free param- 
eters (refer to text for a more detailed description). Here 
A2 = 0.3. The small inset represents A\ vs t/to (squares), 
obtained keeping A2 = 0.3 fixed. The solid line is a power- 
law fit that gives exponent £ = —0.22 ~ —2/9. 



outside the core, where the density reduces drastically and 
the multi-stream flow approaches the single stream limit. 
The squares in the small inset represent, in double log- 
arithmic scale, the amplitude A\ resulting from the fit 
(keeping now A2 = 0.3) , for different times. The power- 
law scaling is clear and the exponent is £ = —0.22 ~ —2/9. 
We will come back on this important point in the end of 
the Section. 

Now, let us consider the system ( ^7|) . Following the dis- 
cussion in Section [l], we perform the change v — > u/ab and 
define the new time variable b. Then the second equation 
of (B7|) reads: 



d b v + (vd x )v = T B (x,v) , 
where T B (x, v) given by: 
T B (x,v) = -J E - p d x [pvar B (x)] = 



9 



■d x [p var Q (x)] 



(41) 



(42) 



ia 4 bp 

Here B indicates that the quantity is expressed in 
Burgers-like variables, i.e. (x,v,b). In fli^ ) the transfor- 
mation w — > ^a(t)v has been considered explicitly. From 
equation pel) we obtain: 



r B (x,v) 



9K (t/t ) 



-d x 



1 



P 



7-1 



a 4 b 2 p 

Inserting the assumption of parallelism 
equation (the third of fl37|)) reads: 

F(t) 



(43) 



the Poisson 



dx P 



AirGa 



d xx u = ~p b bd xx v . 



(44) 



D. Fanelli and E. Aurell: Asymptotic behavior of a stratified perturbation 



9 



Hence equation ( }43| ) is transformed into: 
9Kb(t/to)t Pb 



-d.r 



8 a 4 b 2 P 1+1 
In addition the following relation holds: 



(45) 



1 



1 



Pl +1 



1 



1 — bd T v 



7+1 



and, thus: 



8 u y a 4 - 3 T6 2 Po 



1 



1 — bd r v 



(47) 



where we made use of equation (fL6[). Collecting together, 
equation © reads: 



<9 fc v + (•ue? a .) U = (Jt(t)d x 



1 



1 — 

where pit) is a time dependent coefficient defined by: 



(48) 



PW 8 U 'a^b^Po \toJ 



(49) 



Equation ( f48|) is defined in the inner region where the 
ansatz ( |3^ ) applies and where the gradient of v is negative. 
We stress that, even though obtained in a different man- 
ner, (|48| ) belongs to the same family of equ ations derived 
by ( Buchert, Dominguez & Perez-Mercader 199£), except 
for a slight modification of p{t) and a different interpreta- 
tion of 7. More important, using the constraints imposed 
by the results in section [|, we will provide a precise indi- 
cation of the values of 7 and £. We therefore proceed to 
select one candidate from the whole family d48|). 

Consider equation (48) and apply the following rescal- 
ing for position, velocity and time: 



dl>= I — ) db. (50) 



Equation (|48| ) is transformed into: 

5-2/3 / \A 2 -Ai-<5 



/ \ o—z/ o / \ 



(A)' 



I - 



(*) 



vd x v 

2/3-A!+A 2 



In order to provide a full consistent picture, we require 
the shock interval to shrink in time with the same rate that 
have been shown to hold for the discrete Vlasov equation 
(Q model). Therefore we have to impose: 



a 
2 



2 

3 ' 



(46) 

that implies the following relation between 7 and £: 



I I 

^ = -2^9- 



(54) 



(55) 



Let us now consider the results reported in Fig. |[ 
Assuming £ = —2/9 (see analysis above and small in- 
set in Fig. ||), from equation ( |55| ) one gets 7 = 2/9 that 
is in agreement with the result of the numerical fitting. 
Therefore we are led to conclude that asymptotically, the 
evolution a of multi-stream flow originated by the one- 
dimensional Vlasov-Poisson dynamics, is mimicked by a 
transport equation in the form: 



d b v + vd x v = p(t)d x 



I - 



2/9 



(56) 



where fi(t) cx (t/t )~ 16/9 and b = {t/t ) 2 / 3 . As we will 
show in the next paragraph, the adhesion model can be 
recovered from our ansatz ( |39|) with a special choice of 7. 
Nevertheless both the value of 7 and the consequent rate 
of compression of Ax s hock are not consistent with the nu- 
merical study of the Q model reported in Section [|. We 
are therefore led to conclude that the Burgers' equation, 
in the limit of vanishing viscosity, is only valid as a qual- 
itative test model and that, on the other hand, the more 
quantitative phenomenology is, instead, captured by the 
transport equation (p6|). 

6.1. The adhesion model: limit of validity 

Let us consider the family of equations (|48| ) derived from 
the ansatz (|39|). Assume 7 = — I : this is the only possible 
choice to recover a Burgers like equation, according to our 
approach. In fact equation ([is| ) then reads: 



d b v + (vd x ) v = bp,(t)d xx v 



(57) 



1 ~i 



(51) 



«-4 



By setting Ai = a/2 + 2/3, A 2 = a/2 and 8 = -2/3, 
equation ( pT| ) simplifies: 



drv + -X2V + vd x v — a\d x 



1 



1 — d x v 



(52) 



All the coefficients are now time independent. Hence, 
equation (|5^) develops shocks of constant width, Ax shock ■ 
Transforming back to the old variables, this implies: 



where /i(t) 

There are two major problems with that result. First 
the choice of 7 is not consistent with the simulations 
for the Quintic model. In fact, if 7 < 0, the width of 
the Gaussian profile ( |38| ) decreases moving to region of 
lower density. That is the opposite of what we observed. 
Moreover, since £ = —0.22 ~ —2/9 (see figure ^|), p(t) de- 
cays too fast to agree with the contraction of Aa; detected 
with the discrete Vlasov approach 

Hence, equation (|57|), directly derived from (0), with 



the assumption (39) fails in reproducing the essence of 



A-X s fic 



, x a/2+2/3 / t \ "/ 2+2 / 3 

• Ax shock oc f — J . (53) 



2 Here p(t) oc it/to) 32 / 9 . By a rescaling procedure (anal- 
ogous to the one adopted in the previous section) it can be 



shown that Ax shock oc (t/to) 



-32/9 
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peculiar aspects of the dynamics, that have been shown 
to hold in the Vlasov like picture. Nevertheless, as al- 
ready stressed by (Buchert, Dominguez & Perez-Mercader 



1999), one should note that the viscous like term vanish 
as t — > oo and therefore, at least from a qualitative point 
of view, the Starobinsky conjecture is justified (formally 
replacing v — > in ([fl]) with t — > oo). 

7. Conclusions 

In this paper we discussed the problem of structure 
formation in a three-dimensional expanding Universe, 
focusing on stratified perturbations, in a pressure-less 
medium. This is done by using extensively the Q model, 
a Lagrangian representation that was derived in a recent 
paper (Aurell & Fanelli 2001 ). The Q model is valid in 
the limit where Newtonian mechanics applies and it is 
shown to be equivalent to the Vlasov-Poisson equations 
for N — > oo. 

In particular we investigated the asymptotic behavior 
of an initial smooth perturbation, by measuring the pro- 
gressive contraction of the inner region, compared to the 
overall expansion. Clear power-law scalings are detected 
and a heuristic explanation is provided. Then we derived 
an asymptotic transport equation for the velocity, con- 
sistent with these observations. By means of a combined 
numerical and analytical procedure, we obtained equation 
(pq). Moreover, we showed that the Burgers equation with 
vanishing viscosity, can be directly derived from the ki- 
netic theory, by assuming the ansatz (|39|). Nevertheless, 
equation ([57]) fails in reproducing the correct asymptotic 
scaling observed and, therefore, we are led to conclude 
that the adhesion approach is valid only as an approxi- 
mate model of structure formation. In fact, even though 
Burgers' equation holds outside the shocks, the adhesion 
picture is shown to agree, only at a qualitative level, with 
the correct Vlasov description, in the massive cores. On 
the contrary, inside the shocks, the more quantitative phe- 
nomenology is captured by the transport equation (|5rj|). 
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